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We search for a first-order phase transition gravitational wave signal in 45 pulsars from the 
NANOGravy 12.5 year dataset. We find that the data can be modeled in terms of a strong first 
order phase transition taking place at temperatures below the electroweak scale. However, we 
do not observe any strong preference for a phase-transition interpretation of the signal over the 
standard astrophysical interpretation in terms of supermassive black holes mergers; but we expect 
to gain additional discriminating power with future datasets, improving the signal to noise ratio 
and extending the sensitivity window to lower frequencies. An interesting open question is how well 
gravitational wave observatories could separate such signals. 


Introduction — The search for gravitational waves 
(GWs) spans many orders of magnitude and encapsu- 
lates a plethora of source phenomena. At very-low fre- 
quencies (~ 1 — 100 nHz), pulsar-timing arrays (PTAs; 
[1-3]) aim to detect GWs through the presence of corre- 
lated deviations to radio-pulse arrival times across an en- 
semble of precisely-timed Milky Way millisecond pulsars. 
There are three PTA collaborations that currently have 
decadal-length timing data from an ensemble of pulsars: 
The North American Nanohertz Observatory for Gravi- 
tational Waves (NANOGrav; [4]), the European Pulsar 
Timing Array (EPTA; [5]), and the Parkes Pulsar Tim- 
ing Array (PPTA; [6]). These three, in addition to the 
Indian PTA (InPTA; [7]), are synthesized into the In- 
ternational Pulsar Timing Array (IPTA; Perera et al. 
8). There are also emerging efforts in China (CPTA; 
[9]), as well as some telescope-centered timing programs 
(MeerKAT; [10]; CHIME; [11]). 


The dominant GW signals at such low frequencies fre- 
quencies are expected to be from a cosmic population 
of tightly-bound inspiralling supermassive binary black 
holes (SMBHBs; [12, 13]), producing an aggregate in- 
coherent signal that we search for as a stochastic GW 
background (GWB), and also individual binary signals 
that we attempt to resolve out of this stochastic con- 
fusion background. However, other more speculative 
GW sources in the PTA frequency range include cosmic 
strings |14, 15], a primordial GWB produced by quan- 


tum fluctuations of the gravitational field in the early 
universe, amplified by inflation [16-18], and cosmologi- 
cal phase transitions [19-23], the latter of which is the 
subject this study. 


The most recent PTA results are from NANOGrav’s 
analysis of 12.5 years of precision timing data from 47 
pulsars [24, hereafter NG12], of which 45 exceeded a tim- 
ing baseline of 3 years and were analysed in a search for a 
stochastic GWB [25, hereafter NG12gwb]. NANOGrav 
reported strong evidence for a common-spectrum low- 
frequency stochastic process in its array of 45 analyzed 
pulsars, where ~ 10 of those pulsars are strongly sup- 
portive, most are ambivalent, and a few seem to disfavor 
the process (although not significantly). No evidence for 
the characteristic inter-pulsar correlation signature im- 
parted by GWs was found. At low frequencies the shape 
of the characteristic strain spectrum was well matched 
to a power-law, with an amplitude and slope consistent 
with theoretical models of SMBHB populations. Under 
a model that assumes the origin of the GWB is a popula- 
tion of SMBHBs, the median characteristic strain ampli- 
tude at a frequency of 1/year is 1.92 x 10715. Interpreta- 
tions of this common-spectrum process as a GWB from 
SMBHBs have since appeared in the literature, showing 
that, if it is indeed so, robust evidence of the distinc- 
tive inter-pulsar correlations should accrue within the 
next several years, followed by characterization of the 
strain spectrum and astrophysical probes of the underly- 


ing population [26, 27]. However, the Bayesian posterior 
probability distributions of the strain-spectrum ampli- 
tude and slope are broad enough to entertain a variety 
of different source interpretations, many of which have 
since appeared in the literature [e.g. 28-31]. 

In this Letter we consider gravitational waves produced 
by first-order cosmological phase transitions, both as an 
alternative origin of the common process measured in 
the NANOGrav 12.5 year Dataset [32-39], and as a sub- 
dominant signal to that produced by SMBHBs. The 
frequency range to which NANOGrav is sensitive cor- 
responds to phase transitions at temperatures below the 
electroweak phase transitions of the Standard Model (i.e. 
T < 100GeV). This has led many to consider higher 
frequency GW observatories, such as LISA and LIGO, 
as the dominant instruments to search for phase tran- 
sitions. However, phase transitions may occur at much 
lower temperatures in particular in hidden sectors [40-— 
42]. Hidden sectors/valleys feature rich dynamics, with 
multiple matter fields and forces, independent of the dy- 
namics of the Standard Model. They appear generically 
in top-down constructions like string theory, and in some 
solutions to the so-called hierarchy problem. In many 
cases, they may be difficult to detect via their particle 
interactions with the Standard Model, but gravity is an 
irreducible messenger. In this regard, PTAs provide a 
powerful complementary probe to the dynamics of hid- 
den sectors already being explored through many terres- 
trial, astrophysical and cosmological probes (see Ref. [43] 
for a recent summary). 

Previous studies on cosmological first order phase tran- 
sition in the context of the NANOGrav results were car- 
ried out in [33, 37, 44, 45]. Our analysis presents two 
main novelties compared to these works: first, we prop- 
erly include the relevant noise sources and discuss the 
impact of backgrounds (like the one generated by SMB- 
HBs); second, we discuss how the results are affected by 
the theoretical uncertainties on the GW spectrum pro- 
duced by first order phase transitions. 

The outline of this Letter is as follows. In the next 
section we briefly summarize the signature of GWs from 
the dominant background of SMBH mergers. We then 
dive into the main subject of this Letter, GWs from a 
first-order phase transition, where we discuss the rele- 
vant parameters characterizing the signal. We then carry 
out an analysis with the NANOGrav 12.5 year dataset, 
finding that the data can be modeled in terms of a strong 
phase transition with a transition temperature around 10 
MeV. The dataset and data model for these analyses are 
exactly as described in NG12 and NG12ewb, respectively. 
All common processes (whether interpreted as being of 
SMBHB or phase-transition origin) are modeled within 
the five lowest sampling frequencies of the array time se- 
ries, corresponding to ~ 2.5—12 nHz. Finally, we discuss 
theoretical uncertainties, and compare the PT interpre- 
tation of the data to the standard interpretation in terms 
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TABLE I. Parameters for the gravitational wave spectrum of 
eq. (4). The values of the parameters (a, b,c) in the spectral 
shape of the bubble contribution are reported in Table II. 


of SMBHB finding no strong preference for one over the 
other. 

GW from SMBHBs mergers — Regardless of origin, 
the energy density of GWs as a fraction of closure density 
is related to the GW characteristic strain spectrum by 
[46] 


where Ho is the Hubble constant (set here to be 
67 km/s/Mpc [47]), and the GWB characteristic strain 
spectrum he(f) is often described by a power-law func- 
tion for astrophysical and cosmological sources: 


he(f) = Acws (4) , (2) 


where Agwep is the amplitude at a reference frequency of 
1/year, and a is an exponent that depends on the origin 
of the GWB. For a population of inspiraling SMBHBs, 
this is a = —2/3 [48]. The cross-power spectral density 
of GW-induced timing deviations between two pulsars a 
and b can be written as 


Aawe f a 3 
ab 1272 A yr, (3) 


where y = 3 — 2a = 13/3 for SMBHBs, and Ta» is the 
Hellings-Downs [49] correlation coefficient between pulsar 
a and pulsar b. 

GWs from first-order phase transition — A first-order 
phase transition (PT) occurs when the true minimum of 
a potential is separated from a false minimum by a bar- 
rier through which a field must locally tunnel. This can 
occur in either weakly coupled (where a scalar field tun- 
nels) or strongly coupled (where a vacuum condensate 
corresponds to the scalar field) theory. Such transitions 
are known to proceed through nucleation of bubbles of 
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TABLE II. Comparison of the bubble spectral shape param- 
eters derived using the envelope and thin wall approximation 
[50] (left column), the semi-analytic approach of reference [52] 
(middle column), and lattice simulations [53] (right column). 
For numerical and semi-analytic results the values of the pa- 
rameters depend on the choice of the scalar field potential, 
we report the range of values obtained for the different scalar 
field potentials considered in the above mentioned works. 


true vacuum which, if sufficiently large, expand in the 
background plasma (still in the false vacuum). Collisions 
of these bubbles, as well as interactions between the ex- 
panding bubble walls and the surrounding plasma, can 
be efficient sources of GWs. 

We characterize the phase transition in terms of four 
parameters: 


e T, — the Universe temperature at which the phase tran- 
sition takes place. 

e a, — the strength of the phase transition, defined as 
the ratio of the vacuum and relativistic energy density 
at the time of the phase transition. 

e 2/H, — the bubble nucleation rate in units of the Hub- 
ble rate at the time of the phase transition, H». 

e Vy — the velocity of the bubble walls. 


The three main sources of GWs associated with a first- 
order phase transition are: (i) collisions of bubble walls, 
(ii) collisions of the sound waves generated in the back- 
ground plasma by the bubbles expansion, and (iii) tur- 
bulence in the plasma generated by expansion and col- 
lisions of the sound-wave. However, in this analysis we 
will not include the turbulence contribution as it usu- 
ally is subleading compared to the sound-wave one, and 
also affected by the largest theory uncertanties (see for 
example [54-56] for recent developments). 

The contribution to the total GW spectrum from bub- 
bles and sound waves collisions can be parametrized as 
(54, 57] 


RAUF) =R Alvu) (E=) (F Se) , (4) 


where the prefactor R œ 7.69 x 10-59; '/ accounts for 
the redshift of the GW energy density, S(-) parametrizes 
the spectral shape, and A(v,,) is a normalization factor 
which depends on the bubble wall velocity, vw. The value 


4 


of the peak frequency today, f°, is related to the value 
of the peak frequency at emission, fx, by: 


ns = fa B T gx 1/6 
emare (O EGEE. 
(5) 


where g, denotes the number of relativistic degrees of 
freedom at the time of the phase transition. The values 
of the peak frequency at emission, the spectral shape, 
the normalization factor, and the exponents p and q are 
reported in Table I for all the production mechanisms 
considered in this work. Due to the finite lifetime [58, 59] 
of the sound waves, to derive Qsw eq. (4) needs to be 
multiplied by a suppression factor Y (Tsw) given by [58]: 


T(Tew) = 1 — (1 + 2Tow Ha) 1? (6) 


where the sound-wave lifetime is usually taken to be 
the timescale for the onset of turbulent behaviors in the 
plasma [60]: Tew * R«/Uys, where the average bubble 
separation is given by R, = (87)!/38-'Max(vy, cs) [61], 
and U? X 3kswa/[4(1 + a,)] [60]. 

Generally both the production mechanism contribute 
to the GW spectrum. However, if the bubble walls in- 
teracts with the surrounding plasma most of the energy 
released in the PT is expected to be transferred to the 
plasma so that the sound waves (and possibly the tur- 
bulence) contribution dominates the GW spectrum. An 
exception to this scenario is provided by models in which 
the bubble walls do not interact with the plasma, or 
by models where the energy released in the PT is large 
enough that the friction exerted by the plasma is not 
enough to stop the walls from keep accelerating (run- 
away scenario). However, determining wether or not the 
runaway regime is realized is either model dependent or 
affected by large theoretical uncertainties. Therefore, 
we perform two separate of analyses. A sound-wave- 
only (SWO) analysis, where we assume that the runaway 
regime is not reached and that the sound wave and tur- 
bulence contributions dominate the GW spectrum; there- 
fore we set kọ = 0, and use the results of reference [62] 
to derive Kzy as a function of v,, and a,. A bubble-only 
(BO) analysis, where we assume that the runaway regime 
is reached and that bubble collisions dominate the GW 
spectrum; we then set vy = 1, Ksw = 0 and Kg = 1. 

We conclude this section emphasizing that, despite re- 
cent progress, large theoretical uncertainties still affect 
the prediction of the GW signal produced in cosmologi- 
cal phase transitions. To get an idea of the impact that 
these uncertainties have on our results we will study how 
the BO analysis is impacted by them. Similar, if not 
larger, uncertainties affect the sound wave contribution 
and would impact the results of the SWO analysis. 

Assuming that the stress energy density of the expand- 
ing bubbles is localized in an infinitesimally thin shell 
near the bubble wall (thin shell approximation), and that 
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FIG. 1. In red (blue) the l-o (68% posterior credible level), and 2-0 (95% posterior credible level) contours for the two- 
dimensional posterior distributions in the (T, œ+) plane obtained in the BO (SWO). The BO analysis has been performed 
with the spectral shape computed by using the envelope approximation (left panel), semi-analytic results (central panel), and 
numerical results (right panel). Specifically, we use (a, b,c) = (1, 2.61, 1.5) for the semi-analytic results, and (a, b, c) = (0.7, 2.3, 1) 


for the numerical results. 


it instantaneously decays to zero after two bubbles col- 
lide (envelope approximation), the bubble spectral shape 
can be derived analytically [50, 63]. The spectral shape 
parameters obtained in this way are reported in the left 
column of Table II. To go beyond these approximations, 
3D lattice simulations are needed. These simulations 
are extremely expensive given the hierarchy between the 
large simulation volume needed to include multiple bub- 
bles, and the small lattice spacing needed to resolve the 
thin walls. Because of the relativistic contraction of the 
wall width, this separation of scales becomes increasingly 
large for increasing wall velocities, making it impossible 
to simulate ultra-relativistic walls. However, the GW 
spectrum can be simulated at lower velocities and the re- 
sults extrapolated to larger values. This is the approach 
taken in Refs. [53, 64], where the authors show that at 
high frequencies the GW spectrum is much steeper than 
predicted by the envelope approximation (b ~ 1.4 — 2.3 
depending on the form of the scalar field potential). An 
alternative approach to the problem has been taken by 
the authors of Refs. [52, 65]. In these works a paramet- 
ric form for the evolution of the scalar field during bub- 
ble collisions is found by using two-bubble simulations. 
This parametric form is then used in many-bubble sim- 
ulations to derive the GW spectrum. They also find a 
steeper high frequency slope (b ~ 2.6 — 2.9) compared 
to the prediction of the envelope approximation. Similar 
discrepancies are found at low frequencies, where both 
the numerical and semi-analytic results find a shallower 
spectrum compared to the envelope approximation (see 
Tab. II). To probe the theoretical uncertainty associated 
with the bubble contribution, we will carry out three sep- 
arate BO analyses utilizing each approach and compare 
the constraint on the phase transition temperature and 


strength. 


Results — We now report our results for the BO and 
SWO analyses. For either of them we consider both the 
case where the only GW signal is the one produced by the 
PT, and the one in which the PT signal is superimposed 
to an astrophysical background produced by SMBHB. 
This latter analyses will give an indication of how difficult 
it will be to disentangle a signal from a phase transition 
from the SMBHB background. The prior distributions 
for the model parameters of all these analyses, in addition 
to other noise characterization parameters, are listed in 
Table III. 


The two parameters that we can constrain the most 
are the transition temperature, Tą, and the phase tran- 
sition strength, œ». Their 2D posterior distributions for 
the PT-only searches are shown in Fig. 1. To assess the 
impact of theoretical uncertainties related to the bubble 
spectrum, for the BO analysis we report the results ob- 
tained by using the three different estimates of the bubble 
contribution to the GW spectrum described in the previ- 
ous section (envelope, semi-analytic, and numerical). We 
can see that at the 1-o (68% posterior credible) level all 
the searches prefer a strong PT, a, Z 0.1, with low tran- 


~N 


sition temperature, Tą <10MeV. At 2-0 (95% posterior 
credible) level the posteriors for the semi-analytical and 
numerical results have support at much higher temper- 
atures, while the envelope results still prefer relatively 
low values. The preference for small values of T, at 
the 1-o level can be understood by noticing (see Fig. 2) 
that the data prefer GW spectra that are peaked at fre- 
quencies below the NANOGrav sensitivity window (i.e. 
f? < 10-° Hz). And, by setting 8/H, = 1 in (5), we see 
that this requirement corresponds to Tą X 10 MeV. The 


low-frequency part of the numerical and semi-analytical 


Parameter Description Prior Comments 
White Noise 
Ex EFAC per backend/receiver system Uniform [0, 10] single-pulsar analysis only 
Qk [5] EQUAD per backend/receiver system log-Uniform [—8.5, —5] single-pulsar analysis only 
Ji [s] ECORR per backend/receiver system log-Uniform [—8.5, —5] single-pulsar analysis only 
Red Noise 
Area red-noise power-law amplitude log-Uniform [—20, —11] one parameter per pulsar 
Yred red-noise power-law spectral index Uniform [0,7] one parameter per pulsar 
Phase Transition 
T. [GeV] phase transition temperature log-Uniform [—4, 3] one parameter for PTA 
Ax phase transition strength log-Uniform [—1.3, 1] one parameter for PTA 
H,./8 bubble nucleation rate log-Uniform [—2, 0] one parameter for PTA 
Uw bubble wall velocity log-Uniform [—2, 1] one parameter for PTA 
Supermassive Black Bole Binaries (SMBHB) 
AcwB common process strain amplitude log-Uniform [—18, —14] one parameter for PTA 
YGwB common process power-law spectral index delta function (yews = 13/3) fixed 


TABLE III. Priors distributions for the parameters used in all the analyses in this work. The prior for the bubble wall velocity 
reported in this table is the one used for the SWO analysis, for the BO analyses we use vw = 1 as explained in the text. 


GW spectra is shallow enough that, at the 2-ø level, the 
data can be fitted also by spectra with peak frequen- 
cies above the NANOGrav sensitivity window. The same 
is not true for the envelope results, which have a much 
steeper low-frequency spectrum; this is the reason why 
the 2-o levels of the envelope results deviate substantially 
from the other two. 

In Fig. 2 we show the GWB spectrum predicted for 
the maximum likelihood parameters of PT-only searches. 
To better illustrate our results, and how the different 
parameters and theoretical uncertainties affect the GWB 
spectrum, we release an interactive version of Fig. 2 at 
this link. 

To understand how the inclusion of the SMBHB back- 
ground affects our results, in Fig. 3 we show the poste- 
rior for the parameters a, and Agwp obtained in the 
PT+SMBHB search. As expected, with the inclusion of 
the SMBHB background, the posteriors for a, stretch 
to lower values where most of the signal is provided by 
the SMBHB contribution.! The Bayesian Information 
Criterion (BIC) [66], defined to be BIC = kInn — 2In£ 
where n = 5 is the number of data points in the fre- 
quency space, k is the number of parameters in the model 
and £ is the maximum likelihood, is also computed. For 
the BO searches, the differences in BIC between the 
PT+SMBHB and SMBHB only searches are found to 
be —0.92, 3.04 and 1.89 for the envelope, semi-analytic 


1 The posterior stops around ax œ 0.05 because of choice of the 
prior, otherwise it woiuld extend down to ax ~ 0. 


and numerical results respectively; similarly the BIC dif- 
ferences between the PT-only and SMBHB-only searches 
are —1.82, —3.18, —1.28. For the SWO analysis, the dif- 
ference in BIC between the PT+SMBHB and SMBHB 
only searches is —4.56, while we get —2.19 for the dif- 
ference between the PT-only and SMBHB-only searches. 
We can then conclude that that the PT+SMBHB and 
PT-only models were neither strongly favored nor disfa- 
vored compared to the SMBHB only model [67]. 

A complete set of posteriors for the parameters of 
the PT-only searches (derived by using the semi-analytic 
spectrum for the BO analysis) are shown in Fig. 4. As 
noted previously, at 1-ø level the data prefer a strongly 
first-order phase transition (a, 2 0.1) taking place at 
temperature T, < 10 MeV; while no strong constraints 
on Vw or H.,/ß is observed. We can also notice that the 
higher values of T, allowed in the 2-0 region are accompa- 
nied by slower nucleation rates (large H,./3). We should 
caution, however, that numerical simulation have been 
performed for phase transition strengths up to a, ~ 0.5 
[68], and that our results for a, Z 0.5 are derived by 
extrapolating the results of these simulations. A similar 
remark should be made for H,./3, numerical simulations 
with values of this parameter close to unity have not been 
performed yet. 

Given the low value of Tą, and the strong constraints 
on new physics at such low scales, we expect the phase 
transition to take place in a dark sector with only fee- 
ble interactions with the Standard Model (SM). In order 
to be consistent with the Hubble parameter constraints 
during the era of Big Bang Nucleosynthesis (BBN) [69], 
the energy of this dark sector must be transferred to the 
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FIG. 2. Maximum likelihood GWB fractional energy-density 
spectrum for the BO (red) and SWO (blue) analyses com- 
pared with the marginalized posterior for the free power 
spectrum (independent per-frequency characterization; red 
violin plot) derived in NG12gwb. For the BO analysis we 
show the results derived by using the envelope (solid line), 
semi-analytic (dashed), and numerical (dot-dashed) spectral 
shapes. For the BO analyses the values of (œx, T4) for these 
maximum likelihood spectra are (0.28,0.7 MeV) for the enve- 
lope results, (1.2,3.4 MeV) for the semi-analytic results, and 
(0.13, 14.1 MeV) for the numerical results. While for the SO 
analysis we get (6.0,0.32 MeV). 


SM before the onset of BBN at T ~ 1 MeV. This leaves 
an allowed range of values for the transition temperature 
given by Tą ~ 1 MeV — 100 GeV. The next data release, 
which adds multiple years of observations and extends 
the the sensitivity window to lower frequency, should be- 
gin to resolve the peak of the spectrum or additionally 
shrink the range of allowed values for T.. 


Conclusions — We performed a search for a stochas- 
tic gravitational wave background from first-order phase 
transitions in the 12.5 year NANOGrav dataset. While 
previous NANOGrav analysis found no evidence yet for 
the inter-pulsar correlation signature of a GWB, the ev- 
idence for a common-spectrum process was significant. 
We found that the data can be modeled by a strong 
(a, > 0.1) phase transition taking place at temperatures 
below the electroweak scale. However, the data do not 
show any strong preference between an SMBHB and a 
PT generated signal, but we expect to gain additional 
discriminating power with future datasets, improving the 
signal to noise ratio and extending the sensitivity window 
to lower frequencies. In particular, data from the Inter- 
national Pulsar Timing Array will allow the baseline of 
observations to be significantly extended, and the num- 
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FIG. 3. 1-o (68% posterior credible level), and 2-0 (95% 


posterior credible level) contours for the parameters AcwB 
and a, in the PT+SMBHB search. In red (blue) the results 
for the BO (SWO) analyses. In this figure we have used the 
semi-analytic results for the bubble spectrum. The posteriors 
do not extend to lower values of a, because of our choice for 
the a. prior: log-Uniform [—1.3, 1]. 


ber of monitored pulsars to be greatly expanded. The 
present quality of the data is such that our results are 
not strongly affected by theoretical uncertainties on the 
GW spectral shape. However, methodological improve- 
ments on determining the origin of the GWB spectrum 
will be needed for future datasets in order to separate 
the signal from a first-order PT from the SMBHB back- 
ground, as well as to constrain the microscopic origins of 
the PT. 
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